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this  paper  we  describe^high-order  accurate  Godunov-type  schemes  for 

the  computation  of  weak  solutions  of  hyperbolic  conservation  laws  that  are 

is 

essentially  non-oscillatory .  Her-sho*  that  the  problem  of  designing  such 


schemes  reduces  to  a  problem  in  approximation  of  functions,  namely  that  of 
reconstructing  a  piecewise  smooth  function  from  its  given  cell  averages  to 
high  order  accuracy  and  without  introducing  large  spurious  oscillations.  To 


solve  this  reconstruction  problem  we  introduce  a  new  interpolation  technique 
that  when  applied  to  piecewise  smooth  data  gives  high-order  accuracy  wherever 


discontinuities 


the  function  is  smooth  but  avoids  having  a  Gibbs-phenomenon  at 
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ON  HIGH-ORDER  ACCURATE  INTERPOLATION  FOR 
NON -OSCILLATORY  SHOCK  CAPTURING  SCHEMES 

Ami  Harten* 

1.  Introduction. 

In  this  paper  we  present  an  Interpolation  technique  that  gives  rise  to 
new  high-order  accurate  non-oscillatory  schemes  for  the  numerical  solution  of 
scalar  conservation  lawss 

(1.1a)  u  +  f(u)  =  u  +  a(u)u  “0  ,  -«><x<<*>,  t>0 

w  X  v  X 

(1.1b)  u(x,0)  -  uQ(x)  ,  -»  <  x  <  “  . 

We  assume  the  initial  data  Uq(x)  to  be  piecewise-smooth  functions  that 
are  either  periodic  or  of  compact  support,  and  denote  the  evolution  operator 
of  the  unique  entropy  solution  by  E(t). 

Let  vh(x  ,t)  denote  a  numerical  approximation  to  a  weak  solution  u(x,t ) 
of  (1.1),  such  that  v^  -  v^Xj,^),  Xj  -  jh,  -  nt,  satisfies  the 
conservation  form 

(1.2a)  v^  -  Vj  “  x(^j+1/2  -  *j-1/2^  =  [Eh(T)*v  * 

Here  E^t)  is  the  numerical  solution  operator,  X  -  T/h  and  the 

numerical  flux,  is  a  function  of  2k  variables 

< 1  ’ 2b )  * j+1/2  "  *(vj-k+l'”’,vj+k) 

which  is  consistent  with  (1.1a)  in  the  sense  that 

A 

(1.2c)  f(u,u,...,u)  ■  f(u) 

* 
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It  is  well  known  that  if  the  total  variation  of  the  numerical  solution 

(1.2)  is  uniformly  bounded  in  h  for  0  <  t  <  T 

(1.3)  TV(vh(«,t))  <  C*TV(u0)  , 

then  any  refinement  sequence  h  0,  t  =  0(h),  has  a  subsequence  hj  ♦  0  so 
that 

L1 

(1.4)  vh^ - *  u 

where  u  is  a  weak  solution  of  (1.1).  Furthermore,  if  all  limit  solutions 

(1.4)  satisfy  an  extropy  condition  that  implies  uniqueness  of  the  initial 
value  problem  (1.1),  then  the  numerical  scheme  is  convergent  (see  [4]). 

Our  goal  in  designing  numerical  schemes  is  to  obtain  a  computer  code  that 
is  both  reliable  and  efficient.  To  achieve  reliability  we  would  like  to  use 
schemes  that  are  total-variation-stable  in  the  sense  of  (1.3);  to  get 
efficiency  we  need  high-order  accuracy.  Unfortunately  it  became  clear  in  the 
development  of  shock-capturing  schemes  that  it  is  not  easy  to  satisfy  both 
requirements  at  the  same  time:  Endowing  a  scheme  with  a  property  that  implies 

| 

automatic  control  over  the  growth  of  the  total  variation  of  the  numerical 
solution  may  very  well  lead  to  a  restriction  on  the  order  of  accuracy  of  the 
scheme. 

The  first  successful  attempt  to  achieve  nonlinear  stability  was  to 
require  positivity  of  the  numerical  solution  operator,  which  led  to  the 

1 

development  of  monotone  schemes.  However  the  requirement  of  positivity  auto¬ 
matically  implies  first  order  accuracy  of  the  scheme  (see  [2]). 

The  next  step  in  the  development  was  to  consider  the  larger  class  of 

I 

total-variation-diminishing  (TVD)  schemes,  where  the  numerical  solution 
operator  is  required  to  diminish  the  total  variation  of  any  BV  function  v 

(1.5)  TV(Eh(T)»v)  <  TV(v)  ;  , 

these  schemes  trivially  satisfy  (1.3)  with  C  “  1  (see  [3]).  We  were  able  to 


V  V-V-\'-V.V-V-V  \  VV  v  v-v-v  v-v 


r 


wwwMBrowtxvfezrei vus  --.fa*  r.-?  i ;  * iW _  w„  "v— _  g* 


■B’v  •:» 


■  T7  r  T7T*  V'1 17  ."%  -?  H 


construct  TVD  schemes  that  in  the  sense  of  local  truncation  error  are  high- 
order  accurate  everywhere  except  at  local  extrema  where  they  necessarily 
degenerate  into  first  order  accuracy,  (see  [9],  [3],  [1],  17],  [8]).  The 
perpetual  damping  of  local  extrema  determines  the  cumulative  global  error  in 
the  I^-norm  to  be  of  (1  +  — )-th  order,  i.e.  only  first-order  in  the  maximum 
norm  but  second -order  in  Lj« 

Recently  ([5])  we  went  one  step  further  and  introduced  a  larger  class  of 
non-osclllatory  schemes,  in  which  the  application  of  the  numerical  solution 
operator  to  any  mesh  function  v  is  required  not  to  increase  the  number  of 
local  extrema  (note  that  this  statement  does  not  depend  on  h  nor  on  the 
smoothness  of  v).  Unlike  TVD  schemes,  which  are  a  subset  of  this  class,  non- 
oscillatory  schemes  are  not  required  to  damp  the  values  of  each  local  extremum 
at  every  single  time-step,  but  are  allowed  to  occasionally  accentuate  a  local 
extremum.  Because  of  this  last  property  we  were  able  to  construct  non- 
oscillatory  schemes  that  are  second-order  accurate  also  in  the  maximum  norm. 

In  the  present  paper  we  relax  the  control  over  the  possible  growth  of  the 
total  variation  of  the  numerical  solution  even  further  and  consequently  are 
able  to  design  schemes  that  are  accurate  to  any  finite  order  r.  These 
schemes  satisfy 

(1.6)  TV(E  (t )  *u)  <  TV(u)  +  0(hr+1),  u  e  (J 

n 

where  U  is  the  set  of  functions  that  are  C  except  at  a  finite  number  of 
points  (in  any  finite  interval)  where  they  may  have  a  discontinuity  in  the 
function  or  some  derivative  (  describes  the  generic  form  of  the  "computable" 
solutions  of  (1.1)).  We  shall  refer  to  this  class  of  schemes  as  essentially 
non-osclllatory  schemes.  The  inequality  (1.6)  ensure-  that  the  scheme  does 
not  have  a  Gibbs-like  phenomenon  of  oscillations  that  are  proportional  to  the 
size  of  the  jump  at  a  discontinuity;  the  permissible  increase  in  total 


variation  is  solely  due  to  the  smooth  part  of  the  function.  Unlike  the 
second-order  accurate  non-oscillatory  schemes  of  [5]  which  are  a  subset  of 
this  class,  essentially  non-oscillatory  schemes  may  occasionally  increase  the 
number  of  local  extrema.  However  these  extra  oscillations  are  on  the  level  of 
truncation  error  in  the  smooth  part  and  therefore  can  be  regarded  as  "non- 
essential  oscillations". 

In  the  following  we  shall  present  high-order  Godunov-type  schemes  that 
satisfy  (1.6),  and  show  that  the  design  problem  reduces  to  solving  a  problem 
in  interpolation;  the  latter  is  the  main  topic  of  this  paper. 
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2.  High-order  accurate  Godunov-type  schemes. 

Let 

(2.1)  u(x)  -  ^  2  u(x+y)dy  =  (^‘uXx) 

denote  the  sliding  average  of  u(x).  The  sliding  average  u(x,t)  of  a  weak 
solution  of  (1.1)  satisfies 

(2.2a)  l^-utx/t)  +  ^  (f  (u(x+  j,t) )  -  f(u(x-  ^,t))]  «0  . 

Integrating  (2.2a)  from  t  to  t  +  T  we  get 

(2.2b)  u(x,t+T)  ■  u(x,t)  -  A  (f  ( x+  -|,t/u)  -  f(x-  -^,t;u)] 

where  A  ■  t/h  and 

(2.2c)  f  (x,t;u)  *  jj-  /q  f  (u(x,t+n)  )dn  . 

Thus  the  exact  weak  solution  of  (1.1)  satisfies  the  following  relation  on  the 
computational  mesh 

(2.3)  Uj  -  uj  -  A[f(xj+i/2'Vu)  "  f^xj-1/2'tn,u)l  1 

where  u^n  *  u(Xj,tn)  is  the  cell  average  of  the  solution  at  time  tn. 

Next  we  compare  (2.3)  to  the  numerical  scheme  in  conservation  for m  (1.2) 

A 

where  we  set  vjj  =  u^  .  We  see  that  if  the  resulting  numerical  flux  fj+1/2 
satisfies 

(2.4.)  'V'/2  '  T  t(u(Xj+V2'V")>a’'  +  0(hr> 

then 

(2.4b)  vf1  -  Ujn+1  +  0(hr+1)  , 

provided  that  the  coefficient  in  the  0(hr)  term  in  (2.4a)  is  sufficiently 
smooth.  Relation  (2.4b)  shows  that  in  the  sense  of  cell-averages  the 
truncation  error  of  the  scheme  is  0(hr  )»  i.e. 

(2.5)  u(t+x)  -  Eh(T).u(t)  -  0(hr+1)  . 

We  shall  refer  to  a  scheme  that  satisfies  (2.4)-(2.5)  as  an  r-th  order 
Godunov-type  scheme.  (Note  the  difference  from  Lax-Wendrof f-type  schemes  that 
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are  derived  by  approximating  a  Taylor  expansion  of  the  solution  and  where  the 
truncation  error  is  made  small  in  a  pointwise  sense . ) 

We  observe  that  although  (2.3)  is  a  relation  between  the  cell  averages  of 
the  solution  at  tn  and  tn+1,  the  evaluation  of  the  numerical  flux  in 
(2.4a)  involves  point  values  of  the  solution.  Since 
(2.6)  u(x)  -  u( x )  -=  0(h2) 

we  have  to  devise  a  technique  to  recover  point  values  from  given  cell  averages 
to  a  desired  accuracy,  in  order  to  obtain  Godunov-type  schemes  that  are  more 
than  second-order  accurate.  The  rest  of  this  paper  is  devoted  to  describing 
an  algorithm  for  the  solution  of  this  problem  in  approximation:  Given  {u^} , 
cell  averages  of  u  e  U  (i.e.  piecewise  smooth  with  a  finite  number  of 
discontinuities)  find  R(x;u)  that  reconstructs  u(x)  to  any  finite  order  r 
(2.7a)  R(x;u)  =  u(x)  +  0(hr) 

wherever  u(x)  is  smooth,  and  such  that 
(2.7b)  R(Xjju)  *  uj  • 

(2.7C)  TV(R(mu))  <  TV(u)  +  0(hr)  , 

here  R  denotes  the  sliding  average  of  R. 

It  is  easy  to  see  that  once  we  solve  this  approximation  problem,  the 
Godunov-type  scheme 

(2.8)  vn+1  =  E^tW11  =  A^EdJ-Rf’jv11)  , 

where  A^  is  the  cell-averaging  operator  (2.1)  and  E(t)  is  the  evolution 
operator  of  (1.1),  is  r-th  order  accurate  and  essentially  non-oscillatory  in 
the  sense  of  (1.6). 

To  see  that  the  scheme  (2.8)  is  essentially  non-oscillatory  we  observe 
that  both  Ah  and  E(t)  are  positive  operators  and  therefore  also  TVD. 
Consequently 

(2.9)  TV(Eh(T)*u)  =  TV(Ah»E(l)»R( »;u) )  <  TV(R(»;u)) 
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and  (1.6)  follows  immediately  from  (2.7c). 

Next  we  show  that  (2.8)  is  r-th  order  accurate  in  the  sense  of  (2.4). 

Let  us  denote 

(2.10a)  V  (*,t)  =  E(t)*R(* jvn )  . 

n 

Using  (2.7b)  and  the  fact  that  vn(x,t)  is  an  exact  solution  of  the 
conservation  law  (1.1a)  we  can  rewrite  (2.8)  in  the  form  (2.3),  i.e. 

( 2 . 1 0b )  v.  vj  X^j+1/2  "  ^  j- 1/2^ 

where 

(2.10c)  ^ j+1/2  “T^O  f^Vn^Xj+1/2ft^dt 

(2.4a)  then  follows  from  (2.7a)  and  the  stability  of  entropy  solutions  of 

(1.1). 

Remarks :  (1)  The  scheme  (2.8)  is  the  abstract  form  of  Godunov-type  schemes. 

In  practice  we  use  an  approximation  to  (2.10)  which  is  obtained  by  using  an 
appropriate  numerical  quadrature  for  the  integral  in  (2.10c),  and  using  an 
approximate  solution  operator  to  evaluate  vn(xj+l/2ft)  in  f2*10®)*  For  more 
details  see  [5]  and  [6]. 

(2)  We  initialize  the  computation  by  taking  cell-averages  of  the  given  initial 
data,  i.e. 

(2.11a)  v°  -  Uj0  -  1  f*02  u0(x^)dy  . 

When  we  apply  the  scheme (2. 8)  N  times,  N»t  =  t,  and  the  initial  data 
are  such  that  the  solution  is  smooth,  we  expect  the  truncation  error  (2.5)  to 
accumulate  linearly.  Thus  at  the  end  of  the  time  loop  we  get  (since  t  = 

0(h) ) 

(2.11b)  v^  =  u(Xj,t)  +  0(hr ) 

Since  in  general  we  are  interested  in  pointwise  values  we  output  R(XjjvN) 
which  gives  us  the  pointwise  data  to  the  desired  accuracy: 


(3.1a)  H^Xjju)  =  u(Xj) 

The  interpolant  Hj^xju)  is  a  piecewise  polynomial  function  of  x,  i.e. 

(3.1b)  H^xju)  “  %[f  j+i/2^x,u)  for  xj  <  x  <  xj+l 

where  q^  j  +  1/2  ^ s  a  polynomial  in  x  of  degree  m.  Wherever  u(x)  is 

smooth  we  get  that 

Jc  Jc 

(3.2)  H  (xju)  =  ~  u( x )  +  0(hm+1"k)  ,  0  4  k  <  m  . 

dxk  m  dxk 

(We  use  here  the  standard  convention  that  k  =  0  corresponds  to  the  function 
itself . ) 

The  new  feature  of  this  interpolation  technique  is  that,  although  u  may 
be  discontinuous,  the  interpolant  H^xju)  is  essentially  non-oscillatory  in 
the  sense  that 

(3.3)  TV(H  (mu))  <  TV(u)  +  0(hm+1 )  . 

nt 

This  interpolation  technique  will  be  used  in  the  following  sections  to  develop 
am  algorithm  for  the  solution  of  the  reconstruction  problem  (2.7). 

To  accomplish  (3.2)  we  take  q^  j+1^2(x;u)  to  be  the  m-th  degree 
polynomial  that  interpolates  u(x)  at  the  m  +  1  successive  points  {x^}, 
im(j)  <  i  <  im(j)  +  m,  that  include  Xj  and  Xj  +  -|*  i.e. 

(3.4a)  <lm,j+1/2(xi;u)  =  u(xi)  '  im<  j *  <  1  <  +  m  ' 

(3.4b)  1  -  m  <  im(j)  -  j  <  0  . 

Clearly  there  are  exactly  m  such  polynomials  corresponding  to  the  m 
different  choices  of  i^j)  subject  to  (3.4b). 

We  have  assumed  that  u(x)  has  a  finite  number  of  discontinuities. 

Hence  for  h  sufficiently  small  there  are  at  least  m  +  1  points  of  smooth- 
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ness  between  any  two  successive  discontinuities.  Consequently  if  (Xj,Xj+^) 
is  an  interval  in  which  u(x)  is  smooth,  there  is  at  least  one  choice  of 
im(j)  such  that  u(x)  is  smooth  in  xim(j)  <  x  <  x^ 

Next  we  show  that  any  algorithm  to  assign  im(j)  (3.4)  to  (xj,Xj+^) 
that  has  the  property: 

(3.5)  smoothness  in  (Xj,Xj+1)  ==>  smoothness  in  (x^  (j)»xi  (j)+m^  • 

m  m 

yields  h^xju)  that  satisfies  both  (3.2)  and  (3.3). 

Let  x0  be  a  point  that  has  a  neighborhood  in  which  u  is  smooth. 

Hence  for  h  sufficiently  small  there  is  an  interval  (Xj,Xj+1)  that 
includes  xq  in  which  u  is  smooth.  It  follows  then  from  (3.5),  (3.4)  and 
the  fact  that  qm  j+-j/2  an  interP°^-atin9  polynomial  that  uses  data  from  an 
interval  in  which  u  is  smooth  that 


(3.6) 


dx 


ir  )+,/2tx’a) 


r  u(x)  +  0(hm+1  k)  for  0  <  k  <  m 
dx 


and  Xj  <  x  <  xj+^  l 


this  implies  (3.2). 

We  turn  now  to  show  (3.3).  First  let  us  consider  an  interval  (Xj,Xj+-|) 
in  which  u  is  smooth.  It  follows  immediately  from  (3.6)  with  k  *  0  that 
(3.7) 


TV ,  .  (q  .  ,  ._)  <  TVr  ,  (u)  +  0(hm+1) 

*  j'  j+V  J+l/2 


Moreover  in  the  intervals  in  which  u  is  smooth  and  ~  is  bounded  away  from 
zero,  it  follows  from  (3.6)  with  k  =  1  that  for  h  sufficiently  small 
<5x  Sn  j+1/2  ^  0  such  an  interval  and  consequently 


(3.8) 


TV [>=  ^ »x j+1 1  (qm,  j  +  1/25  TV(Xj,Xj+1] 
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Next  let  us  consider  an  interval  (Xj,Xj+^)  that  contains  a  single 
discontinuity  of  u.  It  turns  out  that  for  h  sufficiently  small  ]+]/2' 
for  any  choice  of  im(j)  in  (3.4),  is  monotone  in  (Xj,Xj+^);  this  implies 
(3.8).  To  simplify  the  argument  let  us  consider  the  case  that  u  is  a  step- 
function  with  the  discontinuity  located  in  (Xj,Xj+^).  In  this  case  u(x^)  = 
u(xi+i)  for  all  i  4  j  and  therefore 

(3‘9a)  %i,j+1/2(Vu)  =  <3m,j+1/2(xi+1,u) 


for  all  indices  i  such  that 


(3.9b) 


i  4  j  and  0  <  i-im(j>  <  m-1  . 


From  (3.9a)  it  follows  that  —  ]+]/2  ^as  a  root  in  each  of  the  m-1 

intervals  (xi/xi+i)  with  indices  i  satisfying  (3.9b).  Since  qffl  ]+\/2 
is  a  polynomial  of  degree  m-1  it  cannot  have  an  additional  root  in 
(Xj,Xj+1)  and  therefore  q^  j+ 1/2  is  strictly  monotone  there. 

We  conclude  from  the  above  analysis  that  an  increase  in  total  variation 
is  possible  only  in  those  intervals  in  which  u  has  a  smooth  local  extremum. 
Since  there  is  only  a  finite  number  of  such  intervals  and  since  the  increase 
there  (3.7)  is  on  the  level  of  interpolation  error,  (3.3)  follows  (see  [6]  for 
a  more  detailed  analysis). 

In  the  following  we  present  an  algorithm  to  assign  i^j)  to  (Xj,Xj+^). 
In  order  to  satisfy  (3.5)  this  algorithm  makes  use  of  the  information  about 
smoothness  contained  in  a  table  of  divided  differences  of  u.  The  latter  can 


be  defined  recursively  by 
(3. 10a) 

(3.10b)  u[xif  . . .  ,xi+Jcl  *  (u (3 


ulXjJ  «  u(x^) 


'xi+k]  ‘  ulxi‘ 


,xi+k-15 )/(xi+k“xi)  * 


It  is  well  known  that  if  u  is  C,,,  in  txi»xi+kl  then 

1  dk 

(3.11a)  ulxi,...,xi+k]  =__£u(!-ifk)  ,  Xi  <  5lrk  <  xi+k  • 
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However  if  u  has  a  jump  discontinuity  in  the  p-th  derivative  in  this 
interval,  0  <  p  <  k,  then 


(3.11b) 


u(xi#...,xi+k)  -  0(h"k+P[^- u])  , 


dx 


rdF  i 

here  [ -  uj  denotes  the  jump  in  the  p-th  derivative. 

d x? 


Our  algorithm  is  recursive:  It  arrives  *.t  Hj.fxju),  which  amounts  to 
assigning  ir(j),  by  successively  evaluating  i^Cj),  k  —  1,...,r.  We  start 
by  setting 

(3.12a)  i^j)  -  j  , 

i.e.  q1  j+1/2  the  first  degree  polynomial  interpolating  u  at  Xj  and 
Xj+1.  Let  us  assume  that  we  have  already  defined  ik(j),  i.e.  qk  j+f/2  is 
the  k-th  degree  polynomial  interpolating  u  at 
(3.13) 


xik(j)'***'xik( j)+k  * 

We  consider  now  as  candidates  for  <Ik+1  j+i/2  the  two  ( k+ 1 ) ~th  degree 
interpolating  polynomials  obtained  by  adding  to  (3.13)  the  neighboring  point 
to  the  left  or  the  one  to  the  right!  this  corresponds  to  setting  ik+i(j)  = 
ik(j)  -  1  or  ik+-|(j)  “  ik(j),  respectively.  We  choose  the  one  that  gives  a 
(k+1)-th  order  divided  difference  that  is  smaller  in  absolute  value: 


(3.12b)  ik+1(j)  -  \ 


ifc(j)  _  1  if  lu[xik( j)-1'*,,,xik( j)+k] I 

<  ^U[Xi.  ( j)' ***'X1.  (j)+k+1]  I  . 


kV*1 


otherwise 


Using  Newton’s  form  of  interpolation  it  is  easy  to  see  that 


V*> 


+k 


(3,14)  qk+1,j+1/2  "  qk,j+1/2  “  u[x 


W» . X.i'i1***11  ‘  i-IJj)  <’t"’‘il 


O'*.’"'  v‘rW.,ri^  •  •  <■.  f.  •*.  f.  T.  Vi  V.  •  *  .  •  .  •  v**  <*  ■  •  -  V  V*  V  ; 

v«v*  v  v\»  >  V  V  v  v  •>  v  V  v  •  •  v  *>  v  v  v  v  *.v  •*.  *•.  •*.  - .  • .  *  .  -  . - 


12 


Since  the  product  in  the  RHS  of  (3.14)  is  the  same  for  both  choices  in 
(3.12b),  we  see  that  our  algorithm  selects  as  q^+.|  j  +  1/2  (k+1)-th 

polynomial  that  deviates  the  least  from  j+1/2  *n  (Xj,Xj+i).  Clearly 
if  h  is  sufficiently  small  it  follows  from  (3.11)  that  the  algorithm  (3.12) 
satisfies  the  requirement  (3.5),  and  consequently  has  the  desired  properties 
(3.2)  and  (3.3).  However  if  h  is  not  small,  (3.14)  shows  that  the  algorithm 
attempts  to  find  the  "least  oscillatory"  polynomial  (subject  to  the  restricted 
choice  in  (3.12b)),  and  thus  is  meaningful  also  in  this  case. 

In  the  next  two  sections  we  describe  two  different  techniques  to  solve 
the  reconstruction  problem  (2.7)  in  terms  of  interpolation.  The  importance  of 
using  the  particular  interpolation  described  in  this  section  is  that  its  non- 
oscillatory  nature  goes  over  to  the  reconstruction  algorithm.  The  non- 
oscillatory  nature  of  R(x;u)  (2.7c)  is  demonstrated  in  section  7  by 
numerical  examples  (analysis  is  presented  in  [6]). 
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Reconstruction  via  the  primitive  function. 

In  this  section  we  apply  a  technique  frequently  used  in  area-preserving 
approximations1  in  order  to  solve  the  reconstruction  problem  in  terms  of 
interpolation . 


Given  cell-averages  u^  of  a  piecewise-smooth  function  u  e  U 


(4.1) 


-  1  fj+1/2  .  .. 

uj  ■  ^  Jxj-1/2 u(y,dl'  '  hj " 


j  ’  Xj+1/2  ‘  Xj-1/2 


we  can  immediately  evaluate  the  point-values  of  the  primitive  function  U(x) 


(4.2) 


U(x)  *  /X  u(y) dy 
X0 


(4.3) 


Since 


(4.4) 


l  K  -  u( 

i=i0  j  j 


xj+ 1/2^  * 


»<*>  *  h D(x> 


we  can  apply  interpolation  to  the  point  values  of  the  primitive  function  (4.3) 
and  then  obtain  an  approximation  to  u(x)  by  defining 
(4.5)  R(xju)  =  Hr(x;U)  . 

We  note  that  this  procedure  does  not  require  uniformity  of  the  mesh. 

The  primitive  function  U(x)  is  smoother  than  u(x)  (by  one  extra 
derivative)  and  therefore  U  e  U.  Hence  we  get  from  (3.2)  that  wherever  u(x) 


is  smooth 


(4.6) 


C  dk  r+1-k 

—  H  (x;U)  -  — —  U(x)  +  0(h  )  ,  0  <  k  <  r  . 

k  r  -  k 

(  dx 


Using  (4.4)  and  (4.5)  we  can  rewrite  (4.6)  as 


(4.7) 


t  l 

$—r  R(xju)  =  ^—7  u(x)  +  0(hr  *)  , 

.  I  ,  l 

dx  dx 


I  thank  Hira  Dyn  from  Tel-Aviv  University  for  pointing  this  out.  A  similar 
approach  was  taken  by  Woodward  [10]  and  Zalesak  [11]. 


<• 


> 

V 

*» 

> 

1 


i 


which  implies  (2.7a)  for  l  *  0. 

We  turn  now  to  study  R(x;u),  the  sliding  average  of  (4.5): 

(4.8)  R(xju)  -  £  J-h/2  R(x+y)<Jy  "  F  [Hr(x  +  U)  “  Hr<X  "  U)1  * 

Denoting  the  interpolation  error  by 

(4.9)  e(x)  -  H^xjU)  -  U(x) 
we  can  rewrite  (4.8)  as 

"R(xju)  ”  ^  [U(x+  ^)  -  U(x-  ^)]  +  [e(x+  ^)  -  e(x-  ^)] 

(4.10)  =  £  u(y)dy  +  IT  le(x+  1]  “  e(x“  1)] 

■  u(x)  +  jj-  (e(x+  j)  -  e(x-  y)] 

Since 


(4.11a) 


e(xj+V2) 


R(Xjiu)  *  Uj 


(4.11b)  e(x)  -  0(hr+1)  , 

we  get  from  (4.10)  that  (4.5)  satisfies  (2.7b),  i.e. 

(4.12a) 

and  that  wherever  e(x)  is  snvooth 

(4.12b)  R(x;u)  ■  u(x)  +  0(h*"  ’)  . 

It  follows  from  (4.12)  that  R(x;u)  is  a  piecewise-polynomial  interpolation 
of  degree  r  of  u.  Note  that  (4.12b)  is  one  order  more  accurate  than  (4.7) 


r+1 , 


with  i  -  0,  which  is  R(xju)  -  u(x)  +  0(h  ). 


i 


i 

i 

i 

• 

i 

t 

t 

i 

» 

f 
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5.  Reconstruction  via  deconvolution. 

In  this  section  we  describe  another  technique  to  reconstruct  a  piecewise 


smooth  function  u  e  U  from  its  given  cell-averages  Uj .  We  assume  that  the 
mesh  is  uniform  and  consider  the  cell-averages  u^  to  be  point  values  of 


u(x). 


(5.1a) 


u<x>  =  iJ-  I X/2  u(x**)dy 


(5.1b) 


Uj  =  u(Xj) 


The  function  u(x)  (5.1a)  was  referred  to  earlier  as  the  sliding  average 
of  u.  Here  we  consider  it  to  be  the  convolution  of  u(x)  with  the 


characteristic  function  of  a  cell  ^ (x) 


(5.2a) 


Vx> 


1/h  for  | x I  <  h/2 

0  for  lx!  >  h/2 


(5.2b) 


u(x)  m  JZ.  “(x-yj^tyjdy  =  (u*<*h)(x) 


In  the  following  we  describe  a  procedure  to  reconstruct  u(x)  up  to 
0(hr)  for  any  finite  rj  this  will  be  referred  to  as  "finite-order 
deconvolution".  Expanding  u(x-y)  in  (5.2b)  around  y  =  0,  we  get 


(5.3a) 


where 


qq  Mr  ^  V  ot> 

-  r  u  ;(x)  (-1)  fh/2  k .  _  v  Uk  <*>.  ^ 

“<*>  -  l  -Ei - —  I-W2  =  1  V"  ,xl 

K*0  K^U 


(5.3b) 


2 

(k+1 ) ! 


k  odd 


k  even 


Multiplying  both  sides  of  (5.3a)  by  h  — -  and  then  truncating  the  expansion 

dx* 

in  the  RHS  at  0(hr)  we  get 


(5.4a) 


-(*),  .  "  v  Uk+I  (k+JU.  .  .  r. 

h  u  (x)  *  I  ah  u  (x)  +  0(h  ) 


r-i-1 


k +1  (k+I), 


Writing  the  relations  (5.4a)  for  X  *  0,...,r-1  in  matrix  form,  we  obtain 


Let  us  denote  the  coefficient  matrix  in  the  RHS  of  (5.4b)  by  C.  This  matrix 
is  upper  triangular  and  diagonally  dominant.  Multiplying  both  sides  of  (5.4b) 
by  C”1  from  the  left  we  get 


Relation  (5.4c)  is  the  essence  of  finite-order  deconvolution. 

Given  u^  ( 5. 1 )  we  apply  the  interpolation  technique  of  section  3  to  u 

t  _  (  £ ) 

and  compute  0(hr)  approximations  to  h  u  (x)  at  Xj  by  taking 
appropriate  derivatives  of  Hm(x/u),  m  >  r-1;  denote  these  approximations  by 
y  0  <  *  <  r-1: 


(5.5a) 


(5.5b) 


D0,J  •  "J 

j  *  h  u^  ^(Xj)  +  0(h  )  ,  X  m  1,...,r-1  . 


Using  D.  .  in  the  RHS  of  (5.4c)  we  obtain  0(hr)  approximations  to 

j 

h*u^(x^)  which  we  denote  by  i*e. 
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(5. 6a) 


)t  j  ■  u(Xj )  +  0(hr) 


(5.6b) 


D.  .  -  h£u(£)(x.)  +  0(hr)  ,  1  <  l  <  r-1  . 

3  J 


Since  C  is  an  upper  triangular  matrix,  D.  .  are  computed  by  back- 

**  j 

substitution:  We  set 


(5.7a) 


Dr-1,j  '  Dr- 1 , j 


and  then  confute  backwards  for  k  =  r-2, ...,0 

r-1 

(5.7b) 


D|t'i  ’  a‘D|l'3  ' 


Finally  we  define 


_  r-1  i  k 

(5.8)  R(x»u)  -  \  —  Dk#j  for  lx_xjl  < 


h/2 


k-0 

It  follows  immediately  from  (5.6)  that 

r-1  1 

(5.9)  R(x»u)  =  l  —  u(k) (x.) (x-x^)k  +  0(hr)  -  u(x)  +  0(hr)  , 

k*0  kl  J  J 

which  implies  (2.7a). 

To  see  that  the  reconstruction  (5.8)  satisfies  property  (2.7b) 
(5.10a) 
we  evaluate 

tr, ..  _  1  rh/2 


R(xyu)  *  u^ 


r(x^ju)  =  ^  J_h/2  R(Xj-y»u)^y 


to  get 


r-1 


r-1 


(5.10b)  R(Xj »u)  -  kt  J*  (  h}  /.h/2  y  dy  ~  °0,j  +  JJ1  °kDk, j  ' 

where  are  (5.3b).  Relation  (5.10a)  then  follows  immediately  from 

comparing  the  RHS  of  (5.10b)  to  (5.7b)  with  k  =  0 

_  r-1 

°0,j  “  °0,j  +  J1  \Dk,j 

and  then  using  (5.5a).  We  note  that  this  result  is  independent  of  the 


particular  values  assigned  to  D.  .  for  1  <  i  <  r-1. 

J 
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Clearly  this  choice  satisfies  (5.5)  for  m  >  r-1. 

(2)  The  finite  order  deconvolution  extends  very  easily  to  the  case  where  the 
mollif ier  fy^y)  is  replaced  by  a  smoother  function  or  a  function  with  a 
different  support. 


6.  The  constant  coefficient  case 


In  this  section  we  consider  the  Godunov-type  scheme  (2.8)  in  the  constant 
coefficient  case 

(6.1a)  U£  +  aUjj  =  0,  a  =  const. 

(6. 1b)  u(x,0 )  =  u0(x) 

The  exact  solution  in  this  case  is  just  pure  translation  with  a  constant 
speed  a 

(6.2)  u(x,t)  -  [E(t)*uQ](x)  =  u0(x-at) 

Consequently  we  can  express  the  numerical  scheme  explicitly  in  the  following 
form: 

(6.3)  Vj+1  =  RtXj-auv11)  ; 
here  r  is  the  sliding  average  of  R. 

The  truncation  error  in  the  sense  of  cell-averages  (2.5)  is  then 

(6.4)  u(x-at)  -  R(x-aTju)  , 

which  measures  how  well  R( • ;u)  approximates  u.  We  have  noted  that 
generally  if  R(»;u)  approximates  u  to  0(hr),  then  R(»;u)  approximates 
u  to  Oth1,  ),  which  corresponds  to  an  r-th  order  accurate  scheme.  To  have 
this  gain  of  one  extra  order  of  accuracy  we  need  spacial  smoothness  of  the 
reconstruction  error;  this  is  evident  from  (2.4a)  and  more  directly  from 
(4.10)  -  (4.11).  When  we  do  not  have  this  smoothness  we  find  that  the  scheme 
is  r-th  order  accurate  in  the  L^-norm  but  not  in  the  maximum  norm  (see  section 
7). 

We  note  that  although  the  problem  to  be  solved  (6.1)  is  linear,  the 
numerical  scheme  (6.3)  is  highly  nonlinear.  The  nonlinearity  enters  through 
the  interpolation,  where  the  stencil  (3.12)  is  chosen  differently  at  each 
point  and  each  time  level,  depending  on  local  smoothness  of  the  numerical 
solution.  In  this  respect  (6.3)  is  conceptually  different  from  standard 


finite-difference  schemes  where  the  stencil  is  arbitrarily  predetermined. 

We  would  also  like  to  point  out  that  this  scheme  breaks  aways  from  the 
somewhat  artificial  notion  of  "upstream  differencing":  The  decision  what 
stencil  to  use  in  the  reconstruction  is  made  on  the  basis  of  smoothness 
considerations  and  has  nothing  to  do  with  the  "direction  of  the  wind".  The 
latter  enters  only  when  applying  E(t)  in  (2.8).  The  resulting  stencil  is  a 
combination  of  the  two. 


7.  Numerical  examples. 

In  this  section  we  present  some  numerical  examples  in  order  to  illustrate 
the  nature  of  the  approximations  described  in  this  paper. 

The  first  set  of  examples  deals  with  smooth  data.  In  Figs,  la  to  If  we 
show  approximations  to  u(x)  =  sin  irx,  -  1  <  x  <  1j  these  were  reconstructed 


via  deconvolution  (see  section  5)  from  the  cell-averages  input 


(7.1a) 


—  —  .  sin(nh/2) 

Uj  =  u(xj)  =  TnhTT)  *  sin<*V 


(7.1b)  Xj  =  -1  +  j  •  |  ,  0  <  j  <  N-1  . 

In  this  example  we  took  N  -  6  and  extended  the  data  outside  [-1, 1]  by 
periodicity.  The  cell-averages  input  is  shown  in  Figs,  la  -  If  by  circles. 
Both  the  reconstructed  approximation  R(x;u)  (5.8)  and  u(x)  *  sin  irx  are 
shown  by  solid  lines. 

Let  us  denote  the  polynomial  degree  of  H(x;u)  (5.12)  by  Pj  and  that 
of  R(x;u)  (5.8)  by  PR: 

(7.2)  Pj  =  deg[H(x;u)]  ,  PR  =  deg[R(x;u)] 

Fig.  la  shows  the  reconstruction  associated  with  the  original  first-order 
accurate  Godunov  scheme:  Pj  *  PR  =  0. 

Fig.  1b  shows  the  reconstruction  associated  with  the  "second-order"  TVD 
schemes  [3]:  Pj  *  PR  =  1*  Comparing  it  to  Fig.  la  we  see  that  the  local 
extrema  are  flattened  in  exactly  the  same  way  -  this  demonstrates  the 
degeneracy  to  first-order  accuracy  there. 

Fig.  1c  shows  the  reconstruction  associated  with  the  second-order 
accurate  non-oscillatory  scheme  of  [5]:  Pj  -2,  PR  =  1.  Comparing  it  to  Fig. 
1b  we  observe  a  considerable  improvement  in  the  quality  of  the 
approximation.  We  also  note  that  the  reconstruction  increases  the  total 
variation  of  the  input  data;  however,  this  increase  is  small  -  it  is  of  the 
size  of  the  truncation  error. 


Figs.  Id/  1e  and  If  show  the  reconstruction  corresponding  to  Pj  =  k/ 

PR  =  k-1  for  k  =  3,4  and  5,  respectively.  We  observe  that  unlike  the 
previous  approximations,  the  reconstruction  here  does  not  go  through  the 
circles;  this  is  a  consequence  of  (2.6). 

The  reconstruction  R(x;u)  is  discontinuous  at  xj+^/2'  s*ze 
the  jump  is  proportional  to  the  reconstruction  error.  As  we  increase  the 
accuracy  of  the  reconstruction  going  from  Fig.  la  to  Fig.  If,  we  observe  that 
the  size  of  the  jumps  becomes  smaller  and  smaller;  in  Fig.  If  the  jumps  are 
hardly  noticeable. 

In  Table  1  we  list  the  L^-norm  of  the  interpolation  error 
| H(x;u)  -  u(x)|  and  of  the  reconstruction  error  |r(x;u)  -  u(x)|  for  a 
refinement  sequence 

(7.3)  N  =  4,  8,  16,  32,  64  in  (7.1)  . 

We  turn  now  to  examine  the  performance  of  the  resulting  Godunov-type 
scheme  (2.8)  in  the  constant  coefficient  case 

(7.4)  uj.  +  ux  “  0,  u(x, 0 )  *  sin  irx  ,  -1  <  x  <  1 

We  input  the  initial  data  in  the  form  of  the  cell-averages  (7.1).  Again  we 

assume  periodicity  in  space,  which  implies  a  period  of  2  in  time  as  well. 

In  Table  2  we  list  both  the  L^-  and  the  L j -error  at  t  =  2  of  the 

scheme  (2.10),  (6.3)  where  R(x;u)  is  reconstruction  via  deconvolution  with 

Pj  and  PR  defined  by  (7.2).  The  calculations  in  this  table  were  performed 

for  the  refinement  sequence  (7.3)  with  \  ■  x/h  *  0.8.  In  addition  to  the 

actual  error  we  also  list  the  quantity  r  which  is  the  computational  order  of 

£ 

accuracy.  Assuming  the  error  to  be  exactly  e^  *  const »h  we  evaluate  r 
for  any  two  successive  calculations  in  the  refinement  sequence  by 

(7.5)  r  -  l°g2(eh/eh/2)  • 
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In  Table  3  we  repeat  the  calculations  In  Table  2  for  the  scheme  (2.10), 
(6.3)  where  now  R(x;u)  is  reconstruction  via  the  primitive  function  (4.5). 
Here  PU  denotes  the  polynomial  degree  of  the  interpolation  H(x;U)  of  the 
primitive  function  (4.2),  i.e. 

(7.6)  Pa  -  deg[H(x;U)J  . 

We  observe  from  Table  2  that  the  schemes  with  Pj  =  k  and  PR  =  k-1  are 
(at  least)  k-th  order  accurate  in  both  La  and  .  Hie  TVD  scheme  with 
P_  =  P_  =  1  is  a  little  bit  better  than  first-order  in  L  and  a  little  bit 

i.  R.  co 

worse  than  second-order  in  Lj. 

We  observe  from  Table  3  that  the  schemes  with  the  odd  order  Pg  =  3 
and  P0  =  5  seem  to  be  third  and  fifth-order  accurate,  respectively,  in  both 
the  and  norms. 

However  when  pu  is  an  even  number,  P„  -  2,  4,  6,  this  order  is 
realized  only  in  the  L^-sense.  In  the  L^-norm  it  seems  to  be  one  order  lower 
in  accuracy. 

Comparing  Table  2  and  Table  3  for  the  fine  mesh  N  =  64,  we  find  that 
the  scheme  from  Table  2  with  Pj  =  k  and  the  scheme  with  PQ  *  k+1  from 
Table  3  seem  to  give  comparable  accuracy. 

We  turn  now  to  examine  the  performance  of  the  various  approximations  when 
applied  to  the  discontinuous  function 

-x  sin(3irx^/2)  for  -1  <  x  <  -  y 

(7.7)  u(x)  -  <  |sin(2TTx)|  for  |x|  <  -1  , 

2x  -  1  -  p  sin(3irx)  for  <  x  <  1 
v  o  3 

which  we  extend  periodically  outside  [-1,1]. 


The  circles  in  both  Figs*  2a  and  2b  show  the  given  cell-average  of  (7.7) 
at  N  ■  40  uniformly  spaced  mesh  points.  The  two  solid  lines  in  Figure  2a 
are  u(x)  and  H(x;u)  with  ^  *  6.  The  two  solid  lines  in  Fig.  2b  are 
u(x)  and  R(x;u)  with  PR  -  5.  We  see  in  both  figures  that  the 
approximations  give  good  accuracy  in  the  smooth  part  and  are  essentially  non- 
oscillatory.  This  is  as  to  be  expected  since  N  ■  40  provides  sufficient 
resolution  of  the  problem  in  the  sense  that  there  are  at  least  7  points  of 
smoothness  inbetween  discontinuities. 

In  Figs.  3a  and  3b  we  repeat  the  calculation  in  Fig.  2  for  N  ■  15.  It 
is  interesting  to  note  that  in  spite  of  the  lack  of  resolution  in  this  case 
both  approximations  are  essentially  non-oscillatory. 

To  compare  the  two  reconstruction  techniques  we  repeat  in  Fig.  3c  the 
same  case  described  in  Fig.  3b,  but  now  using  reconstruction  via  primitive 
function  with  P0  “  6.  The  main  difference  is  in  the  rounding  of  local 
extrema  at  discontinuities  present  in  Fig.  3b»  this  is  due  to  the  min  mod 
operation  in  (5. 12b). 

Hext  we  apply  the  Godunov-type  scheme  (2.10),  (6.3)  to  the  solution  of 
Uj.  +  ux  *  0  with  the  initial  data  (7.7).  As  before  we  initialize  the 
computation  by  taking  cell-averages  of  the  initial  data  and  use  X  *  t/h  *0.8 
In  Figs.  4,  5  and  6  we  present  calculations  performed  with  reconstruction  via 
deconvolution  using  Pj  ■  PR  *  1  in  Fig.  4;  Pj  ■  2,  PR  *  1  in  Fig.  5  and 
Pj  ■  4,  PR  ■  3  in  Fig.  6.  The  results  are  presented  at:  (a)  t  ■  2, 

(b)  t  *  4.  The  circles  in  these  figures  show  the  reconstructed  numerical 
approximation  at  the  mesh  points;  the  solid  line  shows  the  exact  solution. 

We  observe  that  in  all  cases  we  get  non-oscillatory  approximations,  and 
that  the  schemes  are  dissipative  in  time.  The  quality  of  the  results  improves 
with  increasing  formal  order  of  accuracy. 
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In  Fig.  7  we  compare  the  two  reconstruction  techniques.  In  Fig.  7b  we 
repeat  the  same  calculation  as  in  Figs.  4  to  6  using  Pj  *  5  and  PR  ■  4.  In 
Fig.  7a  we  use  reconstruction  via  primitive  function  with  Py  -  6  in  (7.6). 

Unlike  the  situation  in  Figs.  3b  -  3c,  here  the  two  schemes  produce  very 
similar  results. 
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Interpolation  Error  Reconstruction  Error 


Table  1 »  Interpolation  and  reconstruction  L^-errors 
for  u(x)  ■  sin  itx. 


N 

Pl-1,PR-1 

Pi-2,Pr-1 

Pl-3,PR=2 

Pi-4,pr«3 

Pi-5,Pr=4 

Pi-6,Pr»5 

4 

1.86X10"1 

1.86x10_1 

7.39x10"2 

7.39x10"2 

3.17x10"2 

3.1 7x 10"2 

8 

6.85X10"2 

2.84X10"2 

7.63X10-3 

3.1 6x 10"3 

9.37X10"4 

3.88x10"4 

16 

1.87X10"2 

3.72X10-3 

5.36X10-4 

1.07x10"4 

1 . 70x10"® 

3.39x10"® 

32 

4.78X10"3 

4.71x10"4 

3.45x10"5 

3.40x10"® 

2. 76x1 0"7 

2.72x10"® 

64 

1.20X10"3 

5.91x10"5 

2.17x10"6 

1.07x10"® 

4.36x10"9 

2.14x10"10 

4 

2.57X10-1 

2.57X10"1 

1.07X10"1 

1.07x10_1 

4.69x10"2 

4.69x10"2 

8 

1.64X10"1 

6.28X10"2 

1.52x10"2 

5.34X10"3 

1.74X10"3 

5.85X10"4 

16 

4.87X10"2 

1.37X10"2 

1.13x10"3 

2.15x1 0”4 

3.28x10"5 

6.22x10"® 

32 

1.27X10"2 

3.27X10”3 

7. 54x10"® 

7.18x10"® 

5.37x10"7 

5.23x10"® 

64 

3.20X10"3 

8.07X10"4 

7.88x10"® 

2.52X10"7 

8.49x10"9 

4. 16x10"10 

1.00 


*0.000 


0.300 


0.000 


Reconstruction  via  Primitive  function. 
r(xju)  with  Py  =  6  vs.  u(x).  N  =  15. 
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